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Abstract 

In this paper we will study the numerical solution of a discontinuous differential system by a Rosenbrock method. 
We will also focus on one-sided approach in the context of Rosenbrock schemes, and we will suggest a technique 
based on the use of continuous extension, in order to locate the event point, with an application to discontinuous 
singularly perturbed systems. 

Keywords: discontinuous differential systems, Rosenbrock methods, continuous extension, event detection, 
one-sided methods, discontinuous singularly perturbed systems. 
2000 MSC: 60H40, 60H07. 



1. Introduction 

In this paper we will study a class of one-step schemes for ODEs, i.e. the class of Rosenbrock methods, in the 
context of ODEs with discontinuous right-hand side, with particular reference to singularly perturbed discontinuous 
ODEs. The class of problems we are dealing with is generally expressed in the form: 

x' = f{x) = \ f f f\ wheaseJfc } 
1 M^) when x £ R 2 y ' 

for t > t and x(i) = v (see [3J, [19], [40], [E]). The state space W 1 is split (locally) into two subspaces i?i and R 2 
by a surface S such that R" = R\ U £ U R 2 . The surface is defined by a scalar event function h : W 1 — > M, so that 
the subspaces R\ and R 2 , and E, are characterized as 

£ = {x e E n | h(x) = 0} , fii = {xe E™| h{x) < 0} , R 2 = {x e E n | h{x) > 0} . (1.2) 

When solving numerically such a discontinuous system, at each integration step the occurrence of a discontinuity 
is checked. What we do in practice, at a general (n + l)-th step, is to check the sign of h(x n ) ■ h(x n +i): if this 
product is greater than zero, then we continue using the same vector field (/i, or f 2 ) as in the n-th step. Otherwise, 
if the product is less than zero, this means that we need to switch to the other vector field, or to a sliding vector 
field; in the event driven approach, that is the approach we will follow, this switching of the vector field requires 
the accurate computation of the event, i.e. the state in which the event function h vanishes; as a matter of fact, 
our coverage will mainly focus on the problems of the event detection. 

Event-driven methods can be applied just when there are finitely many event points. This class of methods is 
widely used (see (U [TTJ [351 [31] [MJ [35] ) ; it seems to be particularly suited since it has been proved (see for example 
|13j . and Section 3 of this paper) that the order of any method (both explicit and implicit, both one-step and 
multi-step) falls to one when discontinuity occurs and the event is not accurately located. Hence, we could say that 
the event location reduces the stiffness of the problem due to the discontinuity. 

Another important issue of this paper is the so-called one-sided approach. As a matter of fact, when we 



numerically solve a discontinuous system of the form ( 1.1 ), we would wish that the vector fields f\ and/or f 2 would 



be defined also beyond E. Neverthless, in certain situations, the vector field is not defined everywhere because 
the model is designed to be applicable only in certain regions of the state space. This particular type of systems, 
sometimes called systems with model singularities, can be found, for example, in |27l 1171 ITS] . 

As simple example (taken from [27] ) we may consider 
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Xy/1 — t when t < 1 , 

when t > 1 ' ^ j 

where one of the vector fields is not defined when t > 1. In general, every differential system whose right-hand side 
involves square roots, logarithms, inverse trigonometric functions, possesses a singularity, where by singularity we 



mean the region of the state space in which the derivative function fi in (1.1 1 is undefined. The events that occur 
in a neighborhood of a model singularity are often referred to as unilateral events, which means that they should 
be detected without allowing the numerical solution to trespass the event itself. 

For this reason, we cannot use implicit methods, since they require evaluating f(x) at possibly singular endpoints. 
This is the rationale of the approach proposed in Section 3. 

We know from Filippov theory (see for example [12]) what can happen when the solution 'hits' the surface S. 
We will focus on the event location: what will happen after the event is not our present concern. 



Finally, we will focus on a particular class of discontinuous ODEs: the discontinuous singularly perturbed 
systems. The presence of a singular perturbation in a discontinuous system arises in a lot of applications (see 

ii^isniMnira). 

Of course these systems are difficult to solve, because of the discontinuity, and, most of all, because of the singular 
perturbation, which introduces a strong stiffness in the problem. The use of Rosenbrock methods is particularly 
convenient for these systems, because of the low computational effort, and of the good stability properties. This 
will be shown in Section 5. 



2. Rosenbrock methods 

As we shall see in the following, Rosenbrock schemes turn out to be very advantageous in the context both of 
one-sided methods and of discontinuous singularly perturbed problems. 

Rosenbrock methods come out from the linearization of diagonally implicit Runge-Kutta methods (see [23]). They 
preserve good stability properties typical of implicit schemes; on the other hand, an s-stage Rosenbrock method 
requires a lower computational effort, since just s linear systems must be solved. [Frequently, as we shall see later, 
we can handle the method in such a way that the s linear systems to be solved have the same coefficient matrix.] 



For an autonomous system, like the one in (1.1), an s-stage Rosenbrock method has the following expression: 



xi = x + Y,t=i hh 
where cxy, 7^ and 6j are the coefficients of the method, and J = f'(x ). 

Each stage of this method requires the solution of s- linear systems with unknowns hi, and with matrix / — r-faJ. 
Most popular Rosenbrock methods (see, for example, [55]) set -fa = 7, for every i = l,...,s. This position is 
computationally advantageous, since all the matrices are equal and we only need one LU-factorization per step. 
Particularly interesting is the use of Rosenbrock methods in a singularly perturbed system with discontinuous right- 
hand side, as we shall see in section 6. 



2.1. An example of order reduction 

The phenomenon of order reduction in discontinuous differential systems, when the event is not accurately 
located, is well-known in the literature. Gear and 0sterby, for instance, deduced the order reduction in Predictor- 
Corrector methods (see [21]). On the other hand, according to the pioneering work of Mannshardt, (see [26]), Lopez 
and Dieci (see [13] ) have accurately computed the global error for explicit Euler method; they have shown that 
there are two contributions to the O(r) term in the global error: the first depends on the fact that Euler's method 
is a first order method, while the second contribution comes directly from the jump and does not depend on the 
order of the method. 

We will give an example of the order reduction of a 2-stages Rosenbrock method, following the approach of [T5] for 
the explicit Euler method. 

What we want to do is to evaluate the leading term of the local truncation error in the discontinuity interval. For 
simplicity of notation, we assume to = 0, but this assumption is not restrictive. Naturally, localizing assumption 



2 



applies, i.e. xq = x(0); also, we call £1 the instant of the event. 

Our 2-stage Rosenbrock methods, applied to the problem (1.1), reads: 

£1 = xq + + where 

(7 - TjJ)k X = Tf(x ), 

(I - TjJ)k 2 = rf(x + h). 



(2.2) 



method (2.2 1 can be written as 



Calling A the inverse of the matrix (I — T7J), and assuming that both Xq and x + rAfi(xo) are in the 

(2.3) 



3 1 

xi =x + -rAf x {x ) + -rAf 1 (x + TAf x {x Q )) - tA 2 /^). 



A Taylor's expansion of the exact solution up to the first order gives: 

x(t) = x Q + Mt&o) + (r - £i)/ 2 (z(£i)) + 0(r 2 ), 
while the local truncation error (l.t.e.) at t\ — r becomes: 

l.t.e. = bfxixo) + (t- ZMxfa)) + 0(t 2 ) - - \k 2 : 
= £i/l(*o) + (r - £i)/ 2 (z(£i)) + 0(r 2 ) - |^r/i(i ) 
- -AtMxq + AtMxo)) + TA 2 h(x a ) + 0(r 2 ). 



(2.4) 



(2.5) 



A first order approximation of our second-order method is sufficient to highlight the proportionality of the local 
truncation error (in the discontinuity interval) with the jump, by means of a positive proportionality factor that is 
less than r. Thus, denoting by p("frJ) the spectral radius of matrix jtJ, if p{^rJ) < 1, then 



A = I + jtJ + 0(t 2 ). 

On the other hand, if v := Afi(xo), we can expand the term fi(xo + tv) as 

/i(i + rv) = h(x ) + tJ„v, 

where J* is the Jacobian of /1 evaluated in some point in Ri. 
By means of (2.6) and (2.7) we get that (2.5) becomes 

l.t.e. = Ui{xo) + (r - £i)/ 2 (*(£i)) + 0(r 2 ) - ^r/x(z ) 

-^/i(^o)+^/i(^o) + 0(r 2 ) 
= £1/1 (*o) + (r - Si)Mxfa)) + 0(r 2 ) - rh(x ) + 0(r 2 ) 
= (t-£i)[/ 2 (z«i))-/i(so)] +0(t 2 ). 

fi(xo) = h (x (£1)) + (z - x (£1)) ^ (a: (6)) + 0(r 2 ), 
z =z(£i)-£i/i(z (6)) + 0(r 2 ), 



(2.6) 



(2.7) 



Finally, since 
and 
we get that 



l.t.e. = (t- [/ a (a (6)) - /1 (1 (6)) ] + 0(r 2 ). 

We have thus shown the order reduction of this second order Rosenbrock method when applied to a discontin- 
uous differential equation like (1.1). 
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2.2. Continuous extension of Rosenbrock methods 

In general, any numerical method for ODEs provides an approximation of the solution at certain mesh points. 
On the other hand, in certain applications (graphics, delay differential equations, initial value problems with driving 
conditions), these discrete values are not enough. We could need a dense output, i.e., a numerical solution defined 
in each point t in the integration interval, [0, T]. The event location, in the context of discontinuous differential 
equations, is one of the cases in which the idea of continuous extension can be advantageous. Naturally, dense 
output formulas can be found in different ways: first of all, trivially, by piecewise linear or cubic interpolants. 
Nevertheless, for Runge-Kutta methods there are more specific manners: perturbed collocated solutions (see [25]). 
and the classical continuous extension proposed in [43,- It is known that, for the classes of Gauss and Radau 
formulas, these collocation methods have almost half the order of the method itself. Some authors (see [IB] and 
[22] ) add some extra stages to achieve an accuracy of 0(t p ), where p is the order of this method. 
Now, in the context of discontinuous differential systems, dense output for Runge-Kutta formulas have been proposed 
both in |15j . and also in [121 123) . In this latter paper, for example, authors are able to find the event point simply 
by seeking the root of a second order continuous extension of the explicit midpoint rule. In this way, further 
evaluations of function / are avoided, and just a second-order polynomial has to be updated. It is noteworth that 
the continuous extensions of Rosenbrock methods have been proposed also in the context of DAEs: in [32] two 
Rosenbrock methods are proposed: a 4-stage Rosenbrock scheme of order 3, and a 3-stage Rosenbrock method of 
order 2. 

Here, we focus on the continuous extension of Rosenbrock methods, whose theory has been investigated -in the 
smooth case- in [30]. The theory of Ostermann, who threads the same path as Zennaro (see [53]) for continuos 
extensions of Runge-Kutta methods, gives a technique for evaluate the solution outside of the mesh, i.e. for 
approximate x(to + cr),Ver e [0, r]. We define 



X{6) =x + ^2b i (6)k i , for < 9 < I, (2.8) 



a continuos extension of Rosenbrock method (2.1), where the functions bi(0) are polynomials and satisfy: 

i) h(0) = 0; 

ii) fti(l) = b t . 



We point out that (2.8 1 only needs already known facts (the stages ki) and can be evaluated cheaply, since just a 
low-order polynomial in 6 must be updated. 

Let us denote by [•] the function integer part; thus we know, by [30] . that ever y Rosenbrock method of order p 
possesses at least one continuous extension of order q = [^— 1, of the form (2.8); moreover the polynomials bi(o~) 
are defined in theorem 1 of [30], and have degree at most q. Nevertheless, Ostermann's theory does not exclude the 
chance of finding a continuous extension that retains the order of the underlying method. 

As a matter of fact, we are going to use a second order continuous extension of a second-order Rosenbrock 
method. The advantage of this choice, with respect to the choice in [52], is that just a 2-stages method is needed 
to get a second-order approximation. 

The Rosenbrock scheme we use has been introduced in 41 and reads 

3 1 

xi(r) = x + -k x + -k 2 , (2.9a) 

(I-77-J)*i=r/(a;o), (2.9b) 
{I-^TJ)k 2 =Tf(x + k 1 )-2Tki, (2.9c) 

for 7 = 1—^. Now, a second-order continuos extension of this method has been proposed in [33] : 

Xx{6) = x + —±—b 1 (6)k 1 + —±—- b 2 (6)k 2 , (2.10a) 

b x {6) = (9 2 + (2 - 6 7 )0, (2.10b) 
b 2 (0) = 2 - 2 7 0. (2.10c) 
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Figure 1: Example of a problem which presents a singularity 



Remark 2.1. Let us assume we are integrating by method When an event has occurred, i.e. when h(xo)h(xi) < 

0, then we need to compute accurately the state vector x(t) such that h{x(f)) = 0. This is simply done by computing 
the root of the scalar function H(t) := h(xi(r)), where x\{t) is the numerical solution of our Rosenbrock method 
in \2. This computation is accomplished by a classical root-finding routine, such as secant or bisection method. 
Every root-finding routine will produce a sequence 7$. Of course, the computation of each term of the sequence 
requires of updating the internal stages k\ and k2, since the internal stages depend on the step size. This can be 



avoided by using the continuous extension (2.10), that has the great advantage (as we shall see in numerical tests, in 
the last section) of preserving the order of the method. Thus, we are going to compute the root of the new function 
H{9) := h(Xi(9)), for < 9 < I, with a great computational saving. 

Finally, we are going to see in the next section that the continuous extension can be very useful also in order to 
give one-sided conditions. 

3. One-sided Rosenbrock methods 

Further to what we said in the introduction about model singularities, we provide an interesting example, 
proposed in |18j . of a planar two-link robotic manipulator with workspace limitations (see Figure [T]). The dynamics 
of this system are described by the following system of ODEs: 

6>i=wi; (3.1a) 
9' 2 =oj 2 , (3.1b) 

for certain functions wi and u> 2 . Now, we could express 9\ and 9 2 as functions of x and y, which denote the position 
of the two arms in the plane: 

9\ = arctan ( — ) — arccos 1 - ; (3.2) 

W ^ 2h^+y2 J 

y-l x sin{6i) 

9 2 = arctan( t^t) ~ "li (3-3) 



y x — liCos{9i) 

where 



x = l 1 cos9 1 + l 2 cos(9 1 + 9 2 ); (3.4) 
y = hsmO! +Z 2 sin(0i +9 2 ). (3.5) 
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Let us assume, for instance, to be in region R\, thus we are integrating vector field f\. We will consider the case 
in which fi cannot be evaluated outside R± U E. In this situation we will consider one-sided Rosenbrock schemes 
that do not require the evaluation of the vector field /i outside R\ U S. This approach has been proposed in [15] in 
the context of explicit Runge-Kutta methods. The idea is very simple: to give sufficient conditions for which the 
internal stages of a Rosenbrock method can be computed even in presence of a model singularity. For example, for 



the two-stages Rosenbrock method (2.9), this means to give conditions for which fi can be evaluated in xq + ki\ 



this will ensure that x\ can be computed. 

We will give sufficient conditions under which one-stage Rosenbrock methods approach the discontinuity from 



one side, and we will make use of the continuous extension (2.10) in order to give one-sided conditions for the 



method (2.9 1 



One-stage Rosenbrock: general case Assuming b\ = 1, one stage Rosenbrock method reads as: 

X\ — Xq + fci, where 
i 1 =r(/- 7 rJ)- 1 /iW. 

As in [12j . we will assume that there is a positive constant S such that 

hZ{x)h{x) > 8 > 0, VxGi?i. 



(3.6) 



(3.7) 



Now, if x(t) is the solution of (1.1 ) in the region then 



jh{x{t)) = hl{x{t))j t x{t) = hl{x{t))h{x{t)), 



and condition (3.7) implies that the function h monotonically increases along a solution trajectory in R\ (close to 



E) until eventually the trajectory hits £ nontangentially. 

We are going to give an analogous condition in the discrete environment of the numerical solution. We assume to 
be in the interval of the event, i.e. xo G Ri and x± G i?2- In this interval, the continuous function H{a) := h{xi{a)) 
changes its sign in [0,r], then at least one r\ G [0,t] exists, such that H(rj) = 0. 

A sufficient condition for r\ G [0,r] to be the only root of the function H is that the straight line segment x\(a) 
intersects S just once: a sufficient condition for this to be true is exactly the analogous of (3.7), that is 



-^h(xi(a)) = hlix^^-x^a) > 0, Vct G [0,r], 



where X\{a) is defined from the (3.6) as 

xi(a) = x + a (I - 7ctJ)~ 1 /i(^o), Vct € [0,r], 

and, obviously, 

-^xi(v) = (I- iaJ)- l h{x ) + a-^-((I - 7 aJ)- 1 )/ 1 (x ), Va G [0,r]. 
Assuming again that p("faJ) < I, for every a G [0,r], we can write 



(3.8) 

(3.9) 
(3.10) 



(7- 7 aJ)- 1 = £( 7( 7j) fe , 



k=0 
oo 



do 



(i-jojr^Y^ko-x-'hJ)*. 



k=l 



Thus, collecting by the same power of a, (3.10) will become 

d 



da 



Xl (a) = A(xo) + 2 7 aJ/ 1 (x ) + 3 7 VV 2 /i(^o) + 



while 



hJ( Xl (a)) = h^xo + oil + joJ + ^o-ij 2 + ....)/i(zo)); 
truncating the power series up to the second power of a, ( |3.I3[ ) becomes 



(3.11a) 
(3.11b) 

(3.12) 
(3.13) 
(3.14) 
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and then, just a Taylor expansion of (3.14) gives 

hl( Xl {a)) « hj(x ) + (<jfj(x ) + ~fa 2 fJ(x Q )J T )h xx (x Q ). 

With this approximation, and by virtue of (3.12), the scalar product in (3.8) can be written as a power series 
expansion in a: 



hj(x 1 (a))—x 1 (a) = h T x fx + v[f 1 T h xx fx+2 1 h T x Jfx]- 



a 2 [3 7 2 /iJ J 2 /i + ^lfJh xx Jh + 1 fJj~ T h xx f 1 ] 



(3.15) 



where all the functions (both matrices and vectors) of the rig ht-ha nd side are evaluat ed a t Xq. 

Let a k be the coefficient of a k in the power series expansion (3.15); hence inequality (3.8) holds if and only if 



^2a k (T k >0, i.e. ^a k a k > -^a k a k . 



k=0 



k=0 



fc=3 



In general this condition is in practice very difficult to verify. What we are going to do, in practice, is to truncate 
the power series up to the third term; thus, we will require just that the the sum of the first three terms of the power 
series in (3.15) is greater than zero for a = r. This condition can be summarized by the following proposition. 

Proposition 3.1. Assuming that all the following function evaluations are performed at Xq, let there exist three 
constants px, p 2 , and p%, all greater than zero, and r > and sufficiently small, such that 

i) h T x h>p x , (3.16) 

U) flh xx f 1 +2 1 hjjf 1 > - P2 , (3.17) 

Hi) 3 7 2 /£ J 2 f x + 2 1 fJh xx Jf l + ^fJj T h xx f x > -p 3 , (3.18) 

iv) Pi-t P2 ~t 2 P3 >0. (3.19) 

Then the function h(x± (a)) is strictly increasing for every a € [0, r] . In particular, there exists a unique r\ £ (0, t), 
such that h(xi(r])) = 0. 

One-stage Rosenbrock: particular case Another approach is possible if we make the foll owing assumption: the 
matrix (/ — r y<rJ)~ 1 is orthogonal, i.e. (/ — 7crJ)~ 1 = (I — "faJ) T . In this case, formula (3.10) will become: 



da 



xi(a) = A(x ) - 2jaJ T f 1 (x Q ), Va e [0,r] 



On the other hand, 



h I(xi{o-)) ~ hj(x ) + hJ x (x Q )(af 1 (x ) - -ycr 2 J T fx(x )), 



hence, dropping the argument of the functions whenever it is X q, the scalar product (3.8) becomes 

h x (xi(o-))-^Xi(a) = hi fx + afjhljx- 
da 

- lo 2 f1hl x ,^fx - 2 1 ahJj T f 1 - 

- 2 1 a 2 f 1 T h xx J T f 1 + 2 1 2 a 3 f 1 T Jh xx J T f 1 . 
Truncating this product to the second power of a, we get 

hj(xx(a))^-xx(a) =hj fx + a (ft h T xx fx - 2 1 hJj T f 1 )- 

- a 2 (2 1 f 1 T h xx J T fx + ifjhl^fx). 



(3.20) 
(3.21) 

(3.22) 



(3.23) 



Finally, sufficient conditions for (3.8) to be satisfied are 
i) hj(x )fx(x ) > Sx, 

") fx ' {xo)h xx (x )fx(x ) - 2jh x r (x )J T fx{xo) > ~Px, 

Hi) 2jf 1 T (x )h xx (x )J T f 1 (x ) + ^fj (x )hJ x (x )J T fx(x Q ) > -p 2 ; 

iv) Sx - rpx - r 2 p 2 > 0. 

Two-stages Rosenbrock: the continuous extension approach In the previous section we have presented a two-stages 
Rosenbrock method (2.9). For this method, there are two possibilities: 
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Figure 2: Internal stage (xq + fei) £ R2 




Figure 3: Internal stage (xq + fei) E Ri 



l.a) h(x + k x ) < 0, and h(xi) > 0, 
l.b) h(x + fa) > 0. 

In case l.b) -see Figure[2j, since h(xo + k\) > 0, we cannot properly compute x\\ in this situation, a step reduction 
is needed, in such a way to find the value a such that h(x + ki(a)) — 0. Then, if X\{a) is above E, we are back to 
case l.a), with step size a, otherwise we continue integrating. 

In case l.a) -see Figure [3} the sufficient condition which guarantees the uniqueness of a such that h(xi(a)) = is 
(3.8 1. Now, the problem is computing -^xi(a). 



In (2.10) we have presented a continuous extension of method (2.9 1. Since this continuous extension has the 



same order of the method, in the discontinuity interval we could confuse the method itself with its continuous 



extension. Thus, we could replace condition (3.8) with its analogous 



±h(X 1 (0)) = hZ(X 1 (9))±X 1 (9)>O, 



ye e [0,1]. 



(3.24) 



The great advantage with respect to (3.8) is in the computation of ^X 1 (6 I ), since the internal stages ki and k 2 do 
not depend neither on a nor on 9. T 
couple of second order polynomials! 



rus the derivative of Xi{9) with respect to 9 is computed just by deriving a 



But we have just argued that condition (3.8) is essentially equivalent to (3.24). Thus we are able to state the 
following theorem. 

Theorem 3.2. Consider the case l.a), and assume that h(xg + ki) < (we know that, with this appro ach, ki does 
not depend on 9). Denote by c the constant 2(1-27) * n ^ e definition of continuous extension (2.10). Moreover, 
assume that there exist constants 8\ > and p\ > 0, and let r > 0, small enough so that, for the continuous 
extension ( 2.10\ ), the following conditions hold: 

1. ^(Jd(fl))c[(2 - 6 7 )fci - 2 7 fc 2 ] > S 1} for all 6»G[0,1]; 

2. /iJ(X 1 (6»))2c(/ci + k 2 ) > -pi, for all 9 e [0, 1]; 

3. S 2 - rp 2 > 0. 

Then, the function h{Xi{9)) is strictly increasing for 9 G [0,1]. In particular, there exists a unique 77 such that 
h(X 1 ( V )) = 0. 
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Proof For every 9 £ [0, 1], 



±h(X 1 (9)) = hl(X 1 (8))^X 1 (6)>0. 



Now, given the continuous extension (2.101 



= c (29 + 2 - 67) fci + c (29 - 2 7 ) k 2 , 



it results that 



h T x (X 1 (8))^X 1 {9) = h'l(X 1 (9))c( (2 - 67) h - 2 7 fc 2 ) 
+ h£(X 1 (8))2cr f (ki + k 2 ). 
Finally, using hypotheses 1., 2. and 3., we get that 

^jX 1 (6)>6 1 -T Pl >0. 

4. Discontinuous singularly perturbed systems 

Discontinuous singularly perturbed systems are very interesting since they combine both the features of singularly 
perturbed differential systems, and the ones of the differential equations with discontinuous right-hand side. 
We are going to study a singularly perturbed system where the discontinuity involves just the derivative of the 
'slow' component, i.e.: 

y =f(y,z); 
tz = g(y,z). 

Here e is a positive small parameter, and the right-hand side of the first equation is defined in the following way: 

t (., „\ _ / z ) when h (y, z ) < , A 9) 

1 [V ' >~ \ h (V, z) when h(y,z) > ' 

where y = y(t) g R s and z — z(t) € M m , and h : R s x M. m —> R is the event function presented in the beginning. 
Here, we remind some basic concepts of singularly perturbed systems, in the smooth case. For a complete coverage, 
see [53], or or also [23]. We know that the solution of this singularly perturbed system can be written as the 
superposition of the outer solution, a smooth function of the independent variable t, that approximates the exact 
solution for values of t well away from the initial instant t = 0, and the initial layer correction, a rapidly decaying 
function of the stretched time which plays an important role just in the initial e-thick boundary layer: outside 
of it is negligible. In symbols: 

y(t,e) =Y(t,e) + er 1 (t,e) 

z(t,e) =Z(t,e) + C(f,e) 1 J 

where Y(t, e) and Z(t, e) are the outer solutions of slow and fast variable, respectively, while r\ (*, e) and £ (*, e) are 
the initial layer corrections of the slow and fast variable, respectively. 



4-.1. Singular perturbation vs DAE 

When dealing with discontinuous singularly perturbed systems, the choice of many authors -mostly when the 
system is linear in both or cither one of the variables- is to consider the reduced system, i.e. the system 



y = f(y,z), 
= g(y,z). 



(4.4a) 
(4.4b) 



obtained by (4.1 1 by setting e = (see, for example, [21 l2~4l, I3T)] !. 

As a matter of fact, this choice makes sense in particular when variable z can be globally expressed as a function 
of y, in equation (4.4b). Indeed, the fact that equation = g(y,z) admits a global isolated solution with respect 
to z of the form z = go(y), is a pretty strong assumption, but it considerably reduces both the dimension and the 
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stiffness of the problem, since by this hypothesis, the DAE (4.4 1 corresponding to the system (4.1 1 is equivalent to 
the so-called reduced- order model: 



V = f(y,9o{y))- 

Under the further assumption that ReSpec^ 2 , the fundamental problems are: 



(4.5) 



1. is this a good approximation? 



2. does the qualitative behaviour of the starting system (4.1| hold, when considering the reduced system (4.5)? 
Answers 

1. First of all we assume that the event is "far" enough from the initial boundary layer, i.e. that the event occurs 
when the transient phase is over. In this way, we can r eason ably neglect the fast decaying term both in the 
slow and in the fast variable, namely, in the notation of (4.3), n (-, e) and £ (-, e), respectively. Moreover, we 
assume that the events are far enough from each other. 

These assumptions guarantee that the fast system is smooth, i.e. the discontinuity surface is not intersected 
by the solution in the initial layer. Of course, by smooth singular perturbation theory, we know that the 



equation (4.5) is just an 0(e) approximation of (4.3). However, we have to notice that, since e is a small 

~ 4 ), very often an 0(e) approximation is sufficiently 



io- 



parameter (for instance, in some real models it is e : 
accurate, in applications. 

For instance, the example proposed in |24j decouples the system in a slow and a fast subsystem, solves them 
separately, and matches the solution in the border of e-thick boundary layer. 

2. The answer, in general, is negative. A nice example can be found in [BJ. Considering the system 







-sign(2z - yt) 




y'2 




-yi - V2 


(4 


ez' 




yi- z 





simple computations show that sliding cannot occur in the system (4.6), because conditions for attractive 
sliding would be y > |e and y < he, which is impossible, since e > 0. 



On the other hand, plugging e = in system (4.6), we observe that the reduced system exhibits attractive 
sliding in all points of the switching manifold. 

Another interesting example is presented in [35]. Here the discontinuous perturbed system 



x ■ 

a/ 



~s\gn[9x + (1 
- x — y 



0)y) 



(4.7) 



depends on a parametere 9 that can be greater or less than zero. It can be shown that the reduced order 
model has a stable equilibrium point in the origin (0, 0), whereas the perturbed system presents, when 6 < 0, 
an exponentially stable periodic orbit around the origin, switching between the two different vector fields F\ 
and F 2 . This system is the object of the numerical tests presented in the last section. 



4-2. Sliding or crossing 

An interesting feature in the treatment of discontinuous singularly perturbed system is the study of conditions 
for sliding or crossing. Assume we are on the switching manifold, i.e. that x G S. From now on, each function 
evaluation will be accomplished in x G S, so we will drop the argument of each function. 

From Fillipov's theory, we know that the only two situations that guarantee the uniqueness of a solution when 
approaching the discontinuity surface are the following ones: crossing and attractive sliding. Crossing simply 
means that the state vector, coming from one of the vector fields (for instance, f\) "hits" the surface £ and crosses 
it instantly. Attractive sliding means that the state vector is forced to move along £ with a yet to be defined vector 
field. For a complete coverage, see [TJ \5\ [T31 EH SO]. Instead, for the definition of the sliding vector field on the 
intersection of surfaces, see [9] ITO ] IT4" 1 I3T] . 

We know that, if n = n(x) — Wh(x), in the notation of system (1.1), 



• crossing occurs if {n T fi)(n T / 2 ) > 0, 

• sliding (both attractive and repulsive) occurs if (n T f\){n T / 2 ) < 0. 
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Now, let us rewrite these conditions in the singularly perturbed case, i.e. in case of system (4.1)-(4.2|. Define 



Fi 



h 

a 



F? 



h 

a 



Thus, condition for sliding, (whether it is attractive or repulsive) is, naturally, 

{n T F 1 )(n T F 2 ) < 0, 

which means, 



1 



1 



hyfi + -h y gj \ h y f 2 + -h y gj < 0, 
and reordering with respect to the power of e _1 , 

(hyfl) (hyfc) + \ ( (hyh) (h z9 ) + (hyfr) (h^) ) + ^ (h^ < 0, 



from which we get, just multiplying for e , the following inequality: 

{hyh) (hyh) e 2 + ( {hyh) (h z g) + (h y f 2 ) (h z g))e + (h z g) 2 < 0. 



(4.8) 



In this way we have just to examine an algebraic inequality in e. But we know in advance that < e << 1. We 
denote by A the coeffient of e 2 , by B the coefficient of e, and C 2 = (h z g) 2 . Naturally, A, B and C are real numbers, 
since they are just sums and products of scalar products. Thus, (4.8) can be rewritten in the following form 

Ae 2 + Be + C 2 < 0. (4.9) 

Now, it becomes clear that, in the limit for e — > 0, the latter inequality is not satisfied. This means, roughly 
speaking, that sliding is less "likely" than crossing, for the singularly perturbed system (4.1 1. In the following, 
sufficient conditions are given for which sliding and crossing occur. 

Proposition 4.1. Sliding occurs if the following conditions are fullfilled: 



A < 0, B 2 - AAC 2 ^ 0, and e > 



-B - y/B 2 +4\A\C 2 



2A 



Proof If A < 0, then B 2 — 4A C 2 > 0; by Descartes' rule of signs, there will be necessarily a positive root and a 



_ -B+y/Bi+4\A\C 2 



-B-y/B 2 +4\A\C 2 



negative one. Thus, inequality (4.9) holds for e < e\ := v 2A , and e > e 2 := ., ( 

to be greater than 0, so the only choice is e > e 2 (since ei < 0). 



Crossing occurs, from Filippov theory, if an analogous of ( |4.8[ ) holds, with the opposite sign, i.e. 

(hyfi) (hyf 2 ) e 2 + ( (h. y h) (h z g) + (h y f 2 ) (h z g))e+ (h z g) 2 > 0, 



which, in the notation of (4.9), becomes 



Ae 2 + Be + C 2 > 0. 



But e has 

(4.10) 
(4.11) 



Proposition 4.2. Crossing occurs, in system (4-1), if one of the following conditions is satisfied: 

1. A > and B 2 - 4AC 2 < 0; 

2. A > and B > 0, assuming also B 2 - AAC 2 > 0. 
Proof 



1. Under conditions in 1, polynomial in (4.11) has no real roots, and assumes values grater than zero for every 
e in R. 



2. Under conditions in 2, from the Descartes' rule of signs, we know that, in the polynomial of inequality (4.11 ), 
if the signs of coefficients do not change, and if B 2 — AAC 2 ^ 0, then the roots of the polynomial, denoted by 



e\ and e 2 , are both negative. Naturally, inequality (4.11 ) holds in the intervals (— oo, ei), and (e 2 , oo). But we 
know in advance that e has to be such that < e << 1, thus inequality (4.11 ) is satisfied. 
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Remark 4.3. Looking at the left hand-side of inequality (4.. 10), we observe that, if the event function h does not 
depend on y (i.e. if h is such that h = h(z)), then the only allowed behaviour of system (4-1) is crossing. 
Conversely, if the switching function h depends just on y, i.e. h = h(y), then necessarily the terms B and C 2 are 
zero. In this case condition for sliding or crossing is driven uniquely by the sign of A. 



Example As an example of the latter remark, let us consider the following system, which is a modified, discontin- 
uous version of an example given in [30] ■ 



yi= z 

V2 = -sign (2/1) yi 
tz = yz-z- ey 1 



(4.12) 



Obviously, switching function h does not depend on z , since h — h{y) — y%. Thus, both B and C are zero, and, 
being A = 



> 0, then we can state that system (4.12) crosses the switching manifold at the event point. 



4-. 3. Numerical issues 
In subsection 



4.1 



we have discussed if it is convenient or not to approximate the system (4.1 1 with its reduced 



order model (4.4 1. From a strictly numerical point of view, it is well known that the great advantage of Ros enbr ock 
methods with respect to implicit Runge-Kutta methods is the linearity. Now, if we integrated system (4.5) by 



a Rosenbrock m etho d, we would lose this advantage (see [13]). For this reason, we have chosen to integrate the 
original system (4.1), and not the corresponding differential algebraic system. 



stepsize t e Global Error Reduction Factor 



1E-3 


1E-2 


2.180627E-4 




0.5E-3 


1E-2 


1.084872E-4 


2.0100 


0.25E-3 


1E-2 


5.426686E-5 


1.9991 


0.125E-3 


1E-2 


2.713762E-5 


1.996 


0.0625E-3 


1E-2 


1.356825E-5 


1.995 



1E-5 


1E-3 


2.202832E-4 




0.5E-5 


1E-3 


1.102479E-4 


1.9980 


0.25E-5 


1E-3 


5.526992E-5 


1.9947 


0.125E-5 


1E-3 


2.765977E-5 


1.9982 


0.0625E-5 


1E-3 


1.382432E-5 


2.0008 



1E-6 


1E-4 


2.202745E-4 




0.51E-6 


1E-4 


1.107174E-4 


1.9895 


0.25E-6 


1E-4 


5.545491E-5 


1.9965 


0.1251E-6 


1E-4 


2.770196E-5 


2.0018 


0.0625E-6 


1E-4 


1.385162E-5 


1.9999 



Table 1: Global Error and Reduction Factor for first order Rosenbrock method BTT3t, for different values of parameter e 



We have considered the system (4.7) with 9 
method 



-0.9. First of all, we have considered one-stage Rosenbrock 



X\ = x Q + ki 
(I-tJ) h = hixo). 



(4.13a) 
(4.13b) 



Now, 



According to the numerical experiments, we observe that the latter method, applied to system (|4.7|), does not 
lose its order, even reducing the value of the parameter e 
applied to the continuous extension X\{9) 



x 



the event is found by the bisection technique 

Once 



6ki, where ki is the same as in (4.13b) and < 9 < 1 



the first event is localized with the requested tolerance (i.e the state vector is on the sliding surface, with a good 
approximation), we have computed the global error in that point, for different values of the step size. Halving the 
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step size, we have also provided the reduction factor, i.e. the ratio between the global error obtained with step size 
r and the one we got using stepsize J. This confirms that scheme (4.13) behaves like a first order method. 

Similar results are provided for the second order Rosenbrock method (2.9 1, together with its continuous exten- 
sion (2.10), used in the context of event location. Table [2] confirms that the continuous extension (2.9) is a second 
order interpolant for the corresponding method. The computational saving of using the continuous extension is 
definetely more evident for the second order method with respect to first order one. 



stepsize r e Global Error Reduction Factor 



1E-3 


1E-2 


7.880118E-5 




0.51E-3 


1E-2 


2.079523E-5 


3.7893 


0.25E-3 


1E-2 


5.760348E-6 


3.6100 


0.125E-3 


1E-2 


1.434942E-6 


4.0143 


0.0625E-3 


1E-2 


3.581187E-7 


4.0068 



1E-5 


1E-3 


9.336405E-7 




0.5E-5 


1E-3 


2.331221E-7 


4.0049 


0.25E-5 


1E-3 


5.825643E-8 


4.0016 


0.125E-5 


1E-3 


1.456912E-8 


3.9986 


0.0625E-5 


1E-3 


3.648825E-9 


3.9928 



1E-5 


1E-4 


8.310706E-5 




0.5E-5 


1E-4 


2.125030E-5 


3.9108 


0.25E-5 


1E-4 


5.678856E-6 


3.7420 


0.125E-5 


1E-4 


1.427827E-6 


3.9772 


0.0625E-5 


1E-4 


3.658980E-7 


3.9022 



Table 2: Global Error and Reduction Factor for second order Rosenbrock d 2 . Q h for different values of e 



5. Conclusion and future work 

In this paper we have studied some issues about the applications of Rosenbrock methods in the context of 
discontinuous differential systems. We focused on conditions for one-sided Rosenbrock methods, and we showed 
the convenience of using the continuous extension in the context of event location. 

We could carry on our research by studying the integration of sliding vector field by means of any implicit and 
semi-implicit Runge-Kutta scheme. 

Of particular interest cuold be also the study of second order differential equations involving discontinuity just in 
the second derivative. These problems arise frequently in impact mechanics. On the other hand, there is a huge 
literature on smooth second-order differential equations, which could be precious in this context (see, for example, 

bed. 
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